//deff ('y=g(x,t)','y=r*x*(1-x/K)')

function y=h(t)
y=K./((((K-N0)./N0).*exp(-r.*t))+1)
endfunction

//on a u0 qui est une matrice. Les ,: c'est pour
//indiquer que l on a un vecteur

function y=f(t,d)
y=r.*d.*(1-(d./K))
endfunction

function u = eulerexplicite(f,u0,t,n)
lt=length(t)
u(1)=u0
for i=2:lt
t1=t(i)
u0=u(i-1)'
t0=t(i-1)
h=(t1-t0)/n
for j=1:n
u0=u0+h*f(t0,u0)
t0=t0+h
end
u(i)=u0'
end 
endfunction

K=3
r=1
N0=100
xbasc()

t=[0:0.01:20];
//t1=[0:0.5:50];
//t2=[0:0.5:50];
n=100
//plot2d( t,h(t),style=2)
//plot2d( t,eulerexplicite(f,N0,t,n),style=6)
//plot2d( t1,eulerexplicite(f,N0,t1,n),style=6)
//plot2d( t2,eulerexplicite(f,N0,t2,n),style=6)


function v = rungekutta(f,v0,t,n)
lt=length(t)
v(1)=v0
for i=2:lt
t1=t(i)
v0=v(i-1)'
t0=t(i-1)
h=(t1-t0)/n
for j=1:n
v1=v0
v2=v0+h/2*f(t0+0.5*h,v1)
v3=v0+h/2*f(t0+0.5*h,v2)
v4=v0+h*f(t0+h,v3)
v0=v0+h/6*(f(t0,v1)+2*[f(t0+0.5*h,v2)+f(t0+0.5*h,v3)]+f(t0+h,v4))
t0=t0+h
end
v(i)=v0'
end 
endfunction

//plot2d( t,rungekutta(f,N0,t,100),style=3)


//essaie pour comparer les différentes fonctions de l'exercice 1.2
//
function [x]=compare(t,k,m)
lt=length(t)
for i=1:lt
k1(i)=k(t(i))
end
x=k1-eulerexplicite(m,N0,t,n)
endfunction


function [x]=compare2(t,k,m)
lt=length(t)
for i=1:lt
k1(i)=k(t(i))
end
x=k1-rungekutta(m,N0,t,100)
endfunction


//plot2d(t,compare(t,h,f),style=2)
//plot2d(t,compare2(t,h,f),style=3)






function y=l(t,u)
y=r.*u.*(1-(u./q))-((u.^2)./(1+(u.^2)))
endfunction

N1=0.1
//r=0.4
//q=10


//xbasc()
//plot2d(t,ode(N1,0,t,l),style=2)
//plot2d(t,ode(2,0,t,l),style=1)
//plot2d(t,ode(5,0,t,l),style=3)
//plot2d(t,ode(10,0,t,l),style=4)
//plot2d(t,ode(3.7,0,t,l),style=5)

//p132 du livre, condition 1 converge vers 0.4633717
// conver 2, vers 3.900565
// conver 3, vers 5.8739895

function y=b(u)
y=r*u*(1-(u/q))-((u^2)/(1+(u^2)))
endfunction


function [x,iter]=Newton(x0,f,jacf,eps,itermax)
x=x0
fx=f(x)
iter=[0:0.5:10];
c=[0:0.5:10];
g=jacf(x)
while norm(fx)>eps & iter<itermax 
Dx=fx/g
x=x-Dx
fx=f(x)
iter=iter+1
g=jacf(x)
end
if norm(fx)>eps & iter>=itermax then
printf("pas de convergence") 
end
endfunction

function g=jacb(u)
g=r-((2*r*u)/q)-((2*u)/(1+u^2)^2)
endfunction


u0=1.9
//sol=Newton(u0,b,jacb,1e-8,500)

function [x]=matr(u0)
lr=length(k)
lq=length(l)
for i=1:lr
for j=1:lq
r=k(i)
q=l(j)
x(i,j)=Newton(u0,b,jacb,1e-2,200)
end
end
endfunction

//sol=matr(u0);

function u=trac(f,t)
n=length(t)
for i=1:n
x=t(i)
u(i)=f(x)
end
endfunction

r=0.4
q=10
a=[0:0.05:10];
deff ('y=sol(x)','y=Newton(x,b,jacb,1e-10,500)')
//plot2d(a,trac(sol,a),style=4)




function e=nbrptfixe(q,r)
m=4*(q^2-(3*(1+(q/r))))
if m>0 then
	e=3 //il y a 3 points fixe
		elseif m==0 then
			e=2 //il y a 1 point fixe
				else m<0
				e=1 //il y a 1 point fixe
end
endfunction
//sol=nbrptfixe(10,0.4)


function dessin(q,r)
linter1=(q(length(q))-q(1))/(length(q)-1)
linter2=(r(length(r))-r(1))/(length(r)-1)
lr=length(r)
lq=length(q)
for j=2:lr
for i=2:lq
solint=nbrptfixe(q(i),r(j))
if solint==1 then
	xset("color",5)
	xfrect(q(i)-(linter1 /2),r(j)+(linter2 /2),linter1,linter2)
	//pour avoir le cadrillage
			xset("color",1)
			xrect(q(i)-(linter1 /2),r(j)+(linter2 /2),linter1,linter2)
	elseif solint==2 then
	 		xset("color",6)
			xfrect(q(i)-(linter1 /2),r(j)+(linter2 /2),linter1,linter2)
			//pour avoir le cadrillage
				xset("color",1)
				xrect(q(i)-(linter1 /2),r(j)+(linter2 /2),linter1,linter2)
		else solint=3
				xset("color",7)
				xfrect(q(i)-(linter1 /2),r(j)+(linter2 /2),linter1,linter2)
				//pour avoir le cadrillage
					xset("color",1)
					xrect(q(i)-(linter1 /2),r(j)+(linter2 /2),linter1,linter2)
end
end
end
//pour avoir les axes
xset("color",1)
plot2d(0,0,rect=[0,0,q(lq)+1,r(lr)+1])
endfunction
//jaune=3 solutions
//mauve=1 solution
//rouge=1 solution

//plot2d([q(i),r(j+1)],[q(i+1),r(j)],style=7)
//[r(j)]
d=[0:0.5:10];
c=[0:0.5:10];
//soleess=(d(length(d))-d(1))/(length(d)-1)
dessin(c,d)
//plot2d(0,0,rect=[-0.1,-3,3,0.1])
//xrect(0,0,1,2)

